Chapter 1 

Two-fluid barotropic models for powder-snow 
avalanche flows 



Yannick Meyapin, Denys Dutykh and Marguerite Gisclon 



Abstract In the present study we discuss several modeling issues of powder-snow 
avalanche flows. We take a two-fluid modeling paradigm. For the sake of simplic- 
ity, we will restrict our attention to barotropic equations. We begin the exposition 
by a compressible model with two velocities for each fluid. However, this model 
may become non-hyperbolic and thus, represents serious challenges for numerical 
methods. To overcome these issues, we derive a single velocity model as a result of 
a relaxation process. This model can be easily shown to be hyperbolic for any rea- 
sonable equation of state. Finally, an incompressible limit of this model is derived. 



1.1 Introduction 

Snow avalanches represent a serious problem for society in mountain regions. 
The avalanche winter of 1999 attracted a lot of attention to this hazardous natural 
phenomenon jUQID. Further development of mountain regions requires an adequate 
level of avalanche safety. Therefore, avalanche protective measures (deflecting and 
catching dams) become increasingly important |9|. During the same winter, several 
avalanches overran avalanche dams, underlining the need for further research in this 
field. Proper design of protecting structures necessitates profound understanding 
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of the snow avalanches flow and of the interaction process with dams and other 
obstacles ||6l[T3l. 

Natural snow avalanches are believed to consist of three different layers: a dense 
core, a fluidised layer and a suspension cloud. Sometimes the surrounding powder 
cloud is absent and we speak about an avalanche in the flowing regime. Obviously, 
transition boundaries between these layers are not sharp and this classification is 
rather conventional. 

The dense core consists of snow particles in persistent frictional contact [8 |. The 
density is of the order of 300 kg/m 3 and the depth of this layer does not exceed 3 
m. The fluidised regime is characterized by particle's mean-free-paths up to several 
particle's diameters. This dynamics at microscopic level explains more fluid-like 
behaviour at large scales. The density of this layer is in the range of 50 - 100 kg/m 3 
and the height is about 3 - 5 m. To model successfully this kind of flows it is crucial 
to know the complex fluid rheology. Finally, these two interior layers can be covered 
by the powder cloud which is a turbulent suspension of snow particles in the air. The 
density ranges from 4 to 20 kg/m 3 and an avalanche in aerosol regime can reach the 
height of 100 m or more [14]. This flow is driven essentially by turbulent advection 
and particles collisions are unimportant. 

In the present study we are concerned with some questions of powder-snow 
avalanche modelling. Since the interface cannot be defined for this type of flows, 
we choose the modelling paradigm of two-phase flows. In this approach the govern- 
ing equations of each phase are spatially averaged to come up with the description 
of the fluid mixture 171 [TBI. 

It is known ifPfl that the front of such an avalanche can develop the speecQ m/ ~ 
100 m/s. For comparison, the speed of sound co in the air is about 300 m/s. It means 
that the local Mach number Ma can reach the value of 

Ma := ^ «0.33. 

Hence, compressible effects may become important. That is why, we begin our ex- 
position with a compressible model. Then, we gradually simplify it to come up with 
an incompressible one at the end of the present article. The goal is achieved by 
taking the limit as the Mach number tends to zero. 

The present article is organized as follows. In Section fTTSl we present a barotropic 
compressible two-phase model with two velocities. Then, this model is simplified in 
Section [T31 using a velocity relaxation process. The incompressible limit of resulting 
system is derived in Section 11.41 Finally, several conclusions and perspectives are 
drawn out in Section [T31 



1 When we estimate the Mach number magnitude, the particle characteristic velocity should be 
taken. However, this information is not easily accessible and we took the maximum front velocity. 
It can lead to some overestimation of the Mach number. 
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1.2 Two-phase flow modelling 

Let us consider a domain Q. C E 3 where a simultaneous flow of two barotropic 
fluids occurs. All quantities related to the heavy and light fluids will be denoted 
by + and — correspondingly. In view of application to snow avalanches, one can 
consider the heavy fluid of being constituted of snow particles and the light fluid is 
the air. When the mixing process is extremely complicated and it is impossible to 
follow the interface between two fluids, the classical modelling procedure consists in 
applying a volume average operator lT7l [T5ll . Thereby, we make appear two additional 
variables a (x,t), x G Q which are called the volume fractions and defined as: 

« x,f := lim i— -1, 
\da\-to \dii\ 
xedn 

the heavy fluid occupies volume d£2 + C d£2 and the light one the volume dQ~ C 
dQ (see Figure fTTt such that 

\dQ\ = \dQ + \ + \dQ-\. (1.1) 

From the relation (II. lb it is obvious that OC + (x,f) + or (x,f) = 1, Vx e £2. 

After performing the averaging process, one obtains two equations of mass and 
momentum conservation: 

5,(a ± p ± ) + V-(a ± p ± u ± ) = 0, (1.2) 
<5 ; (a ± p ± u ± ) + V-(a ± p ± u ± ®u ± ) + a ± V;5 = V- (a ± T ± ) + a ± p ± g, (1.3) 

where p ± (x,f),u ± (x,f ), T ± (x,f) are densities, velocities and viscous stress tensors 
of each fluid respectively. Traditionally, the vector g denotes the gravity accelera- 
tion. We assume that both fluids share the same pressur^l p = p ± (p ± ) and equations 
of state of each phase fulfill minimal thermodynamical requirements: 




Fig. 1.1 An elementary fluid volume dQ occupied by two phases. 



2 In general, this kind of assumptions is reasonable, since relaxation processes will tend to equili- 
brate the system when time evolves. 
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p ±(p±) >0 , dp± d ^ ] > 0, for p ± >0. (1.4) 



In order to obtain a well-posed problem, governing equations ( 11.21 , d 1 - 3b should be 
completed by appropriate initial and boundary conditions. 

If we assume both fluids to be Newtonian, the viscous stress tensor t takes the 
following classical form: 

T ± =A ± trB(u ± )I + 2/i ± B(u ± ), trB^) = V • u*, (1.5) 

where I := {8ij)i<ij<3 is the identity tensor, D(u) := — (Vu^ +'(Vu ± )^ is the 

deformation rate and A ± , /i are viscosity coefficients. For ideal gases, for example, 
these coefficients are related by Stokes relation A ± + |/x ± =0. In application to 
powder-snow avalanches, viscosity coefficients X , jU should be understood in the 
sense of eddy viscosity. 

Remark 1. From physical point of view, presented here model ( 11.2b . ( 11.3b is far from 
being complete. For example, one could supplement it by capillarity effects in the 
Korteweg form. Also we omited all the terms which model mass, momentum and 
energy exchange between two phases. Generally, their form is strongly dependent 
on the physical situation under consideration. 

Remark 2. Since we do not consider the total energy conservation equation, the flu- 
ids are implicitly assumed to be barotropic. In the absence of viscous stresses T , 
the flow is isentropic. This simplification can be adopted provided that important 
energy transfers do not occur. Non-isentropic flows are considered in [fTTl . 

Remark 3. While considering two-phase flows, it is useful to introduce several addi- 
tional quantities which play an important role in the description of such flows. The 
mixture density p and mass fractions m are naturally defined as: 

p(x,r) := a + p + + a~p~ > 0, V(x,r) e Q x [0,T], 

± a± P ± + . - i 
m := m +m =1. 

P 

The total density p is assumed to be strictly positive everywhere in the domain £2 . 
Hence, the void creation is forbidden in our modeling. 

Important quantities p, w ± will appear several times below. 



In principle, one could use equations dl.2b . dl.3b to model various two-phase 
flows. However, this system remains quite expensive for large scale simulations 
required by real-life applications. The major difficulty comes from the advection 
operator associated to model ( 11.2b . ( 11.3b which can be non-hyperbolic |3][15). In 
the next section we will derive a simplified two-fluid model which is proposed as a 
candidate for powder-snow avalanche compressible simulations. 
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1.3 Velocity relaxation 



We would like to reduce the number of variables in the system ( 11.2b . (11.3b . The 
main idea is to introduce the common velocity field for both phases. For this pur- 
pose, we will introduce a relaxation term to the momentum conservation equation 
(O) : 

^(a ± p ± u ± ) + V-(a ± p ± u ± (giu ± ) + a ± Vp = V-(a ± T ± ) + a ± p ± g±-(u + -iO, 

(1.6) 

where K = 0(1) is a constant and e is a small parameter which controls the mag- 
nitude of the relaxation term. Physically this additional term represents the friction 
between two phases. In the following, we are going to take the singular limit as 
the relaxation parameter £ — * 0. This is achieved with Chapman-Enskog type ex- 
pansion. In this way, we constrain velocities u ± (x,f) to tend to the common value 
u(x,f). This technique has been already successfully applied to the Baer-Nunziato 
model &\ in lfl2l 

The first step consists in rewriting the governing equations (11.2b . (11.6b in the 
quasilinear form. To shorten notations, we will also use the material time derivative 
which is classically defined for any smooth scalar function <p(x,f) as 

d ± m d(D , _ 
dt dt 

Lemma 1. Smooth solutions to equations ( 17.21 ), ( 17.61 ) satisfy the following system: 

a±d ^ +P±{c > )2d -af + a± ^) 2V ' U± = °' (1-7) 
a ± p ± c L^_ + a ±v p = v • (« ± t ± ) + a ± p ± g±-(u + - vT), (1.8) 



where (cf) 2 : = |^ 



dt e 

represents the sound speed in each phase ±. 



Proof. This result follows from direct calculations. First of all, we remark that the 
mass conservation equation dl.2l ) can be rewritten using the material derivative as 
follows: 

^^ + aVV-u ± =0. (1.9) 
dt 

Using equations of state p = p ± (p ± ), we can express the density material derivative 
in terms of the pressure and the sound speed: 

d ± p ± 1 d ± p 
dt ~ (cf) 2 dt 

Now, it is straightforward to derive equation ( 11. 7b from ( 11. 9b . 

Finally, if we multiply equation ( |1.9b by and subtract it from the momentum 
conservation equation ( 11.6b . we will get desired result ( 11.8b . 
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Equations ( 11.71 ), ( 11.81 ) can be also recast in the matrix form which is particularly 
useful for further developments: 



A(V £ )^ +B(V £ )VV £ = V • T(V £ ) +S(V £ ) + 



(1.10) 



where we introduced several notations. The vector V £ represents four unknown 
physical variables V £ := '(p,0£ + ,u + ,u~) and 2& := ' (d t p,d t a + ,d t u + ,d t u~) and 
VV e := ' (Vp, Va+, (-V)u+, (• V)u") . Matrices A(V e ) and B(V e ) are defined as 



fa+ p + (c+) 2 
or -p-(cj f 
a+p+I 

\ orp-y 



B(V e ) := 



/a+u+ p+(c+) 2 u+ a+p+(c+) 2 I 
a u -p~(c~) 2 \T orp~(c7) 2 I 
a+I a+p+u+ 

\ a"I a"p _ u _ y 

In these matrix notations the size of zero entries must be chosen to make the multi- 
plication operation possible. 

On the right hand side of ( 11.10b , the work of viscous forces is denoted by symbol 
V-T(Ve) := r (0,0,V-T+,V •<!;-). The source term S(V e ) :='(0,0,a+p + g,a~P~g) 
incorporates the gravity force and R(V e ) := '(0,0, K"(u + — u~), — k(u + — u - )) con- 
tains the relaxation terms. 

Since we expect the limit V £ — > V to be finite as e — > 0, necessary the limiting 
vector V lies in the hypersurface R(V) = 0. In terms of physical variables, it implies 
u + = u . Consequently, we find our solution in the form of the following Chapman- 
Enskog type expansion: 

V e = V + eW + 0(e 2 ). 

After substituting this expansion into ( 11.10b and taking into account that R(V) = 0, 
at the leading order in e one obtains: 



dV 

A(V)— +B(V)VV = V -T(V) + S(V) +R'(V)W, 
at 



(1.11) 



where 



R(V) 



/0 \ 
00 

00 Kl -K-I 
\0 — JCl fCl / 

Henceforth, we make a technical assumption of the presence of both phases in 
any point x G £2 of the flow domain. Mathematically it means that < a + < 1 . 
Since a + + a~ = 1, the same inequality holds for a~. Otherwise, the relaxation 
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process physically does not make sense and we will have some mathematical tech- 
nical difficulties. 

Under the aforementioned assumption, the matrix A(V) is invertible. Hence, we 
can multiply on the left both sides of ( 11.1 lb by PA 1 (V) where the projection matrix 
P is to be specified below: 

p^+PA- 1 (y)B(y)vy = PA- 1 (y)v-T(y)+PR'(y)w+PA- 1 (y)S(y), 

at 



where R'(y ) := A (y )R'(y) and has the following components 



(1.12) 



R'(V) = 






K 



/oo 

00 
00 

—I 



\ 



-I 



V 



a-p- ) 



The vector of physical variables V has four (in ID) components '(/?, a + ,u,u) 
and only three are different. In order to remove the redundant information, we will 
introduce the new vector U defined as U :—' (p,a + ,u). The Jacobian matrix of this 
transformation can be easily computed: 



/l 0\ 
1 
00 1 

Vooiy 



J - du 

In new variables equation ( 11.12b becomes: 

PJ^ + PA" 1 (C7)B(C/)JVC7 = PA~'(t/)V • T(C7) + PR'(U)W + PA 1 (U)S(U). 
at 

(1.13) 

Now we can formulate two conditions to construct the matrix P. First of all, the 
vector W is unknown and we need to remove it from equation (11-1 3b - Hence, we 
require PR'(y) = 0. Then, we would like the governing equations to be explicitly 
resolved with respect to time derivatives. It gives us the second condition PJ = I. 
The existence and effective construction of the matrix P satisfying two aforemen- 
tioned conditions 

PR'(y) = o, pj = i, 

are discussed below. Presented in this section results follow in great lines lfl2l . 

We will consider a slightly more general setting. Let vector V G K" and its re- 
duced counterpart U € R"~ k , k < n. In such geometry, R'(V) G Mat„,„(R), J € 
Mat„ ) „_ J t(M) and, consequently, P e Mat„_^„(M). Here, the notation Mat,„,„(M) de- 
notes the set of m x n matrices with coefficients in M. We have to say also that from 
algebraic point of view, matrices R'(V) and R'(V) are completely equivalent. Thus, 
for simplicity, in the following propositions we will reason in terms of R'(V). 
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Lemma 2. The columns of the Jacobian matrix J form a basis o/ker(R'(V)). 

Proof. If we differentiate the relation R(V) = with respect to U, we will get 
the identity R'(V)J = 0. It implies that range(j) C ker(R'(V)). By direct com- 
putation one verifies that dim range (R'(V)j = k. From the well-known identity 
range(R'(V )) ker(R'(V )) = M", one concludes that dimker(R'(V)) =n — k. But 
in the same time, the rank of J is equal to n — k as well. It proves the result. 



Theorem 1. We suppose that for all V, range(R'(V)) nker(R'(V)) = 
exists a matrix P G Mat„_^„(M) such that PR'(V) = and PJ = 



{0} then it 



Proof. Hypothesis range(R'(V )) n ker(R'(V)) = {0} implies that range(R'(V )) 
ker(R'(V)) = R". From Lemma it follows that range(j) = ker(R'(V)). Thus, 
the space W can be also represented as a direct sum range(R'(V)) 0range(j). 
We will define P to be the projection on ker(R'(V)J = range(j). Since obviously 
R'(V ) G range(R'(V)) and J G range(j), we have two required identities: PJ = I n _£ 
and PR'(V) = 0. 

Now, in order to compute effectively the projection matrix P, we will construct 
an auxiliary matrix D(V ) = [J 1 , . . . ,J"~ k ,I l , . . . ,/*], where J' is the column i of the 
matrix J and {I 1 , . . . ,I k } are vectors which form a basis of range (R'(y)) . We remark 
that PD(y) = 0]. Lemma |2] implies that the matrix D(V) is invertible. Thus, 
the projection P can be computed by inverting D(y): 

p=[i„_,,o]-D- 1 (y). 



Let us apply this general framework to our model (11.121 . where n = 4 and k = 1 . 
The matrix D(y) and its inverse D~' (V) take this form: 



/1 \ 



D(V) = 



1 







K 


00 1 






a+p+ 




K 


00 1 





D- 1 (y) = 



i, 

a-p- J 



(\ 
1 
m+I 

trvmT p 





mT\ 

m + m~ p 
K 



where m are mass fractions defined in Remark[3] 

Now, the projection matrix P can be immediately computed: 

'10 
10 
,0 m+ImT, 



Finally, after computing all matrix products PA~' (f/)B(t/)J, PAr'(f/)V • T(C7), 
PA _1 (f/)S(f/) present in equation ( 11.13b . we obtain the desired single velocity 
model: 
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-^+u-Vj> + pc?V-u = 0, (1.14) 

da + 

+ u- Va + + a+cr<5V-u = 0, (1.15) 

dt 

p^+p(u-V)u = pg + V-T, (1.16) 

where p = a + p + + a~p~ is the mixture density and c\ is the sound velocity in the 
mixture which is determined by this formula: 

2< p+p-(c+) 2 (c7) 2 



pc s := 



a~p + (c+) 2 + a+p-(c s ) 2: 
and 8 is given by 

5;= P + (cf) 2 -p-(c 7 )2 



a-p+(c+) 2 + a+p-(c s ) 2 ' 

Finally, t := A trO(u)I + 2pO(u) is the viscous stress tensor of the mixture. Vis- 
cosity coefficients A, p are naturally defined as 

A := a + A + + a~A~, p := a + p + + a~\iT. 

Equations ( 11.14b - ( 11.161 ) can be recast in the conservative form which is more 
convenient for numerical computations and theoretical analysis. To achieve this pur- 
pose, we replace the pressure p in d 1 . 14b by p ± using the equation of state: 

^+u.V P ± + /4v.u = 0. 

dt (cf)2 

The last equation is multiplied by a , the second equation ( 11.15b is multiplied by 
p ± and we sum them to come up with two mass conservation equations. Trans- 
formation of the momentum conservation equation ( 11.16b is straightforward. The 
resulting conservative system takes this form: 

<9 ( (« ± p ± ) + V-(a ± p ± u) =0, (1.17) 
<9,(pu) + V-(pu®u) + Vp = V-T + pg. (1.18) 

These equations represent a barotropic version of the four-equations model pro- 
posed in J31|5]. 

It can be shown that the advection operator of the model ( 11.17b . ( 11.181 ) is hy- 
perbolic for any reasonable equation of state ( |1.4b . Moreover, this system contains 
fewer variables which allow more efficient computations required in practice. 
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1.4 Incompressible limit 

The main scope of this paper is certainly around compressible two-fluid models. 
However, we decided to derive an incompressible limit of the single velocity model 
( 11.171 i, ( 11.1 8b for the case when acoustic effects should be filtered out. The presence 
of acoustic waves represent, for example, a major restriction for the time step, if an 
explicit scheme is used. 

For the sake of simplicity, we will neglect dissipative effects which do not af- 
fect the acoustic wave propagation. Thus, in this section we consider the following 
system of equations: 

<9,(a ± p ± ) + V-(a ± p ± u) =0, (1.19) 
p<?,u + p(u- V)u + Vp = pg. (1.20) 

(1.21) 

For convenience, we rewrite equation (11.181 1 in nonconservative form. 

In order to estimate the relative importance of various terms, we introduce dimen- 
sionless variables. The characteristic length, time, and velocity scales are denoted 
by £, to and Uq respectively. For example, I may be chosen as the diameter of the 
fluid domain £2, to is the biggest vortex turnover time and Uq is the typical flow ve- 
locity. The density and the sound velocity scales are chosen to be those of the heavy 
fluid, i.e. Pq and Cq s correspondingly. Since we are interested in acoustic effects, 
the natural pressure scale is given by p ( ^(c ( ] i ) 2 . If we summarize these remarks, de- 
pendent and independent dimensionless variables (denoted with primes) are defined 
as: 

•=-?■ '=-5- u: =zv (p ): =i?' p:= ^i?' 

Remark 4. There is nothing to do for the volume fractions a*, since this quantity is 
dimensionless by definition. 

After dropping the tildes, nondimensional system of equation becomes: 

Std,(a ± p ± ) + V-(a ± p ± u) = 0, (1.22) 

Stpd,u + p(u-V)u+ J^V/7 = -^pg, (1.23) 
Ma z Fr z 

where several scaling parameters have appeared: 

e 

• Strouhal number St := . In this study we will assume the Strouhal number 

U t 

£ 

to be equal to one, i.e. ?q = — . 

Uo 

• Mach number Ma := — r- which measures the relative importance of the flow 
speed and the sound speed in the medium. 
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• Froude number Fr := -^2= compares inertia and gravitational force. This param- 

eter will not play an important role in the present study. 

All physical variables a*, p ± , p and u are expanded in formal series in powers of 
the Mach number: 



<p = <po + Ma<pi -f Ma (fa 



(pe{a ± ,p ± , P ,u}. (1.24) 



Formal expansion dl.241 i is then substituted into the system dl.221 i, J 1 .23b . At the 
orders Ma ~ 2 and Ma ~ 1 , we obtain 

Vp = Vpi = 0. 

In other words, po — po(t) and p\ = p\ (r) are only functions of time. At the order 
Ma° we get the following system of equations: 

dt(<4Po) + V ' (^PiN) = 0, 

po<3 r u + po(u ■ V)u + Vn = -3 Pog, 

Fr 



(1.25) 
(1.26) 
(1.27) 



where by n we denote p2- 

Using the same asymptotic expansion dl.24t , one can show that at the leading 
order we keep usual relations between densities and volume fractions: 



a Q+ a Q =1 ) P0 = «0 + P0 + + «0 PQ 



(1.28) 



In order to investigate the behaviour of p^, we will invert the equation of stat43 
p ± = p ± (p) = (jr^ )~ l (p) and expand it in powers of Ma : 



P ± (p)=P ± (po)+Ma 



dp- 



dp 



p\ +Ma 



PO 



.(dp* 
\ dp 



2„± 



P2- 



d 2 P 



On the other hand, from (11.24b we know that 



Po 



2„± 



dp 2 



Po 



pf +0(Ma- 



P =Pc f+ Ma Pj T+Ma L pl 

Matching these expansions at two lowest orders shows that p^ are functions only 
of the time variable: 



Po ± =P ± (Po(0)=:'-o ± (0, Pf = 



dp 



dp 



Po(') 



Pi(t)=:rf(t). 



It is possible to show that p^ l are just constants. Consider the Gibbs relation which 
reads 



3 The function p = p ± (p ± ) is invertible since it is a strictly increasing function > 0. 



dp_ 
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T ± ds ± = de ± - -^dp ± . 

Since we consider isentropic flows, ds^ 1 = and, consequently, the Gibbs relation 
takes a much simpler form: 



It can be shown by considering the total energy conservation equation [11], that 
the internal energy naturally scales with Uq. After dividing ( 11.291 by dt and 
switching to dimensionless variables, equation ( 11.291 ) takes the following form (after 
droping the primes): 

de ± _ p dp ± 
~dT ~ Ma 2 (p±) 2_ rfr' 

Expanding e in the series ( 11.24b and looking at two leading terms, leads to the 
desired result: 

d Pai + 
— — =0 => Pn i = const. 
dt Ul1 

The incompressibility condition V • uo = is obtained by summing up mass conser- 
vation equations (11.251 and taking into account relation (11.28b . 

If we summarize all developments made above and switch back to dimensional 
variables, the resulting incompressible system will become: 

a,a ± + Va ± u = 0, (1.30) 
Vu = 0, (1.31) 
p<?,u + p(u- V)u + V7T = pg + V-T, (1.32) 

where we dropped the index and added again dissipative effects. Viscous stress 
tensor T is still defined by expression ( 11.51 ), as in compressible case. In this case, we 
can speak about two-fluid Navier-Stokes equations. This system of equations ( 11.30b 
- ( 11.321 ) is much easier to solve numerically than its compressible analogue ( 11.171 ), 
( 11.181 ). In particular, this simplification is due to removed stiffness of acoustic waves. 



1.5 Conclusions and perspectives 

In this study we presented several barotropic two-fluid models which can be used 
for numerical simulation of powder-snow avalanche flows. One of the main objec- 
tives of this paper was to reveal the connection between barotropic models with 
single and two velocities. The extension to more general fluids is in progress ifTTl . 

Our exposition began with compressible two-phase model ( 11.21 ), ( 11.31 ) possessing 
two velocity variables. Then, using a relaxation process, we constrained the system 
to have a common velocity for both phases. Mathematically it was achieved with a 
Chapman-Enskog type expansion. Resulting model ( 11.17b , ( 11.181 ) is hyperbolic for 
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any reasonable equation of state fll.4| i. Finally, two-fluid Navier-Stokes equations 
( 11.30b - (11.32t were derived as an incompressible limit of the single velocity model 
(fllTt . CCH . 

Hence, we presented three different two-fluid models which are related by for- 
mal derivation procedures. Simplifications made above, represent a good trade-off 
between accuracy and computational complexity. The final choice should be made 
after determining the flow regime and main goals of the simulation. 

We did not incorporate yet any turbulence modeling. In this study we were fo- 
cused essentially on the advection operators. However, it is obvious that the physical 
flow under consideration is fully turbulent in its aerosol part |[T4ll . As the first phys- 
ical approximation, turbulence effects can be taken into account by adding eddy 
viscosity terms and, thus, by modifying the viscous stress tensor z. It will be done 
in future studies. 



Acknowledgement 

The authors would like to acknowledge the University of Savoie for the PPF 
grant linked to the project: "Mathematiques et avalanches de neige, une rencontre 
possible?". The support from the Research network VOR (Professors Jacky Mazars 
and Denis Jongmans) and Cluster Environnement through the program "Risques 
gravitaires, seismes" is also acknowledged. 

We would like to thank Professor Carmen de Jong for interesting discussions 
around snow avalanches. Special thanks go to our colleagues and friends Didier 
Bresch and Celine Acary-Robert for their continuous help and support. Finally, the 
second author thanks Professors Jean-Michel Ghidaglia and Frederic Dias for intro- 
ducing him to the beautiful field of two-phase flows. 



References 

1. Ancey, C, Bain, V., Bardou, E., Borrel, G., Burnet, R., Jarry, E, Kolbl, O., Meunier, M.: 
Dynamique des avalanches. Presses polytechniques et universitaires romandes (Lausanne, 
Suisse) (2006) 

2. Baer, M., Nunziato, J.: A two-phase mixture theory for the deflagration-to-detonation tran- 
sition (DDT) in reactive granular materials. International journal of multiphase flow 12(6), 
861-889 (1986) 

3. Bresch, D., Desjardins, B., Ghidaglia, J.M., Grenier, E.: On global weak solutions to a generic 
two-fluid model. To appear in Archive for Rational Mechanics and Analysis (2009) 

4. Dias, E, Dutykh, D., Ghidaglia, J.M.: A two-fluid model for violent aerated flows. Submitted 
to Comput. & Fluids (2008) 

5. Dutykh, D.: Mathematical modelling of tsunami waves. Ph.D. thesis, Ecole Normale 
Superieure de Cachan (2007) 

6. Dutykh, D., Acary-Robert, C., Bresch, D.: Numerical simulation of powder-snow avalanche 
interaction with an obstacle. Submitted to Applied Mathematical Modelling (2009) 

7. Ishii, M.: Thermo-Fluid Dynamic Theory of Two-Phase Flow. Eyrolles, Paris (1975) 



14 



Y. Meyapin, D. Dutykh & M. Gisclon 



8. Issler, D.: Experimental information on the dynamics of dry-snow avalanches. In: K. Hutter, 
N. Kirchner (eds.) Dynamic Response of Granular and Porous Materials Under Large and 
Catastrophic Deformations, vol. 11. Springer, Berlin (2003) 

9. Johannesson, T., Gauer, P., Issler, P., Lied, K.: The design of avalanche protection dams. Tech. 
rep., European Commission (2009) 

10. Lied, K.: Satsie: Avalanche studies and model validation in europe. Tech. rep., European 
Commission (2006) 

11. Meyapin, Y., Dutykh, D., Gisclon, M.: Velocity and energy relaxation in two-phase flows. 
Submitted to Eur. J. Mech. B/Fluids (2009) 

12. Murrone, A., Guillard, H.: A five equation reduced model for compressible two phase flow 
problems. J. Comput. Phys. 202, 664-698 (2005) 

13. Naaim-Bouvet, E, Naaim, M., Bacher, M., Heiligenstein, L.: Physical modelling of the inter- 
action between powder avalanches and defence structures. Nat. Hazards Earth Syst. Sci. 2, 
193-202 (2002) 

14. Rastello, M., Hopfinger, E.: Sediment-entraining suspension clouds: a model of powder-snow 
avalanches. J. Fluid. Mech. 509, 181-206 (2004) 

15. Rovarch, I.M.: Solveurs tridimensionnels pour les ecoulements diphasiques avec transferts 
d'energie. Ph.D. thesis, Ecole Normale Superieure de Cachan (2006) 



